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Propagation  of  3-D  Beams  Using  a  Finite-Difference 
Algorithm:  Practical  Considerations 

Alan  H.  Paxton 

Air  Force  Research  Laboratory,  Directed  Energy  Directorate, 

Kirtland  Air  Force  Base,  NM  87117-5776 

Abstract 

We  discuss  practical  aspects  of  the  use  of  the  finite-difference,  alternating-direction  implicit  (ADI),  numerical 
solution  of  the  paraxial  wave  equation,  for  the  free-space  propagation  of  light  beams.  Results  of  calculations 
solving  the  finite-difference  equations  are  compared  with  fresnel-integral  solutions.  Calculations  are  for  round 
beams,  but  the  field  is  represented  in  cartesian  coordinates.  Modes  for  empty  unstable  resonators  are  also 
obtained  using  the  finite-difference  algorithm  and  are  compared  with  fresnel-integral  solutions. 

Keywords:  beam  propagation;  wave  optics;  paraxial  wave  equation;  finite  difference 

1.  INTRODUCTION 

A  finite-difference  method  may  be  useful  for  beam  propagation  in  the  numerical  simulation  of  lasers.  The 
representation  of  the  field  at  many  planes  along  the  resonator  axis,  in  the  gain  medium,  may  be  necessary  due  to 
temperature  variation,  stress-induced  birefringence,  or  dependence  of  the  refractive  index  on  the  light  intensity. 
A  finite-difference  algorithm  can  be  computationally  efficient  for  such  cases.  For  3-D  propagation  in  cartesian 
coordinates,  the  alternating-direction- implicit  method  is  unconditionally  stable  (for  empty-space  propagation). 
The  values  for  the  field  at  the  next  plane  are  obtained  by  solving  a  tri-diagonal  matrix  equation  for  each  row  and 
column  of  the  field  grid.  The  application  of  boundary  conditions  at  the  edges  of  the  x,  y  grid  is  a  consideration, 
where  z  is  the  propagation  direction.  We  use  a  zero-field  boundary  condition  with  apodization  of  the  held  in 
the  vicinity  of  the  boundary.  For  the  symmetrical  test  cases  shown  below,  we  could  have  used  the  boundary 
conditions  for  the  complex  Held  amplitude 

du 
dx 
du 
dy 

Instead,  the  entire  Held  array  was  included  on  the  grid  because  the  code  was  written  to  solve  nonsymmetrical 
problems. 

Cartesian  coordinates  are  often  preferable  to  cylindrical  coordinates,  even  for  treating  round  beams.  All  of 
the  cells  of  the  grid  have  the  same  dimensions,  so  wavefront  aberrations  of  the  same  scale  can  be  treated  at  any 
location  in  the  beam.  Moreover,  round  beams  do  not  usually  have  cylindrical  symmetry  in  high-average-power 
non-fiber  lasers.  The  symmetry  is  broken  due  to  misalignment  of  optical  elements,  flow  in  gas  lasers,  or  the 
cooling  geometries  of  solid-state  lasers. 

Rensch  used  finite-difference  propagation  in  a  simulation  of  a  CO2  laser.1  More  recent  work  on  3-D,  finite- 
difference  optical  propagation,  including  non-paraxial  methods,  was  reviewed  and  augmented  by  Bekker.2 

2.  FINITE  DIFFERENCE  APPROXIMATION  TO  THE  PARAXIAL  WAVE 

EQUATION 

The  paraxial  wave  equation  can  be  written 

V2u  -  2 ik^-  =  0 


along  the  line  x  =  0 

(1) 

along  the  line  y  =  0 

(2) 

(3) 
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The  problem  is  posed  as  an  initial-value  problem,3  with  the  z  spatial  coordinate  replacing  the  time  coordinate  of 
problems  for  which  “initial-value”  is  taken  literally.  Applying  the  ADI  method  to  this  equation,  two  half-steps 
are  applied  successively.  The  first  half-step  is  implicit  in  x  and  the  second  half-step  is  implicit  in  y. 


1(7  \( Ui+1t2  2 Ui+1/2  I  I  (vi  -  2U  4V  'il  7/i+1//2  -I-  7/*  —  f) 

o  [\uj- l,fc  zuj,k  +  Uj+l,k>  +  \Uj,k-l  Zuj,k  +  uj,k+i)i  uj  k  -\-Uj  k-  u 


where 


Az 

2fc(As)2 


and  As  =  Ax  =  Ay 


(5) 


Index  i  corresponds  to  the  z  coordinate,  j  corresponds  to  x  ,  and  k  corresponds  to  y.  The  second  half-step  is 
taken, 


—  la  U  i+l/2  0  i 

~  2uj 


])  +  -  H+k  +  -  utk  +  <41/2  = 0 


More  compact  notation  is  sometimes  used 


(^ui+1/2  +  sy)  -  ui+1/2  +  u*  =  0  (7) 

^Y~(sy+1  +  5lui+1/2)  -  ui+1  +  u?;+1/2  =  0  (8) 

Each  half-step  requires  the  solution  of  a  tri-diagonal  matrix  equation  for  each  row  or  column  of  the  field  array. 
For  the  results  presented  here,  we  set  the  complex  field  amplitude  u  to  zero  on  the  boundary.  The  matrix 
equation  for  the  first  half-step  is, 


(9) 


where  A  =  —1;  B  =  2(i/a  +  1);  Ri  =  Uj  k_1  +  2{i/a  —  1  )ulj  k  +  Uj  k+1,  and  the  transverse  field  grid  has  m  by  m 
points.  For  all  i,  j,  and  k,  we  set:  ulo  k  =  0;  0  =  0;  ulm+l  k  =  0;  u*-  m+1  =  0.  The  accuracy  depends  on  the 

value  of  <j,  with  numerical  error  occurring  for  large  values  of  a. 

The  boundary  condition  would  cause  reflection  of  optical  flux  arriving  at  the  edge  of  the  grid.  To  prevent 
this,  a  smoothed  aperture  was  included.  The  complex  amplitude  was  clipped  by  a  slightly  apodized  circular 
aperture  after  each  half-step  propagation.  The  complex  amplitude  was  multiplied  by  a  cosine  window  function 
with  a  width  of  5  points.  The  aperture  was  located  so  that  the  cosine  window  was  barely  inside  the  grid  at  its 
closest  points.  For  the  calculated  results  presented  here,  the  calculational  grid  was  square,  of  width  twice  the 
initial  diameter  of  the  circular  beam. 


3.  COMPARISON  OF  SOLUTIONS  FOR  PROPAGATION  OF  TOP-HAT  BEAMS 

We  compare  finite-difference  solutions,  for  propagation  of  a  beam,  with  fresnel-integral  solutions.  The  fresnel 
integral  was  solved  using  a  kernel  averaged  algorithm  developed  by  Dwight  Phelps.4,5  A  beam  with  uniform 
intensity  and  planar  wave  fronts  passes  through  an  initial  circular  aperture,  the  radius  of  which  is  1.  The  aperture 
edge  is  usually  slightly  apodized.  The  beam  with  a  slightly  smoothed,  top  hat  distribution  propagates  to  a  final 
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plane  corresponding  to  a  given  fresnel  number.  Abrupt  aperture  edges  would  create  field  variation  for  which  the 
transverse  laplacian  could  not  be  represented  by  the  difference  equations,  so  smoothing  of  the  field  over  several 
points  might  be  required  at  the  initial  beam  aperture.  For  the  propagation  of  an  apodized  top-hat  beam,  unless 
the  smoothing  distance  is  smaller  than  or  comparable  to  the  outer  fresnel-zone  width,  d,  the  apodization  will 
affect  the  solution.  The  width  is  given  by 


d/a  =  1  — 


n  —  1 
n 


1/2 


where  a  is  the  aperture  radius  and  n  is  the  fresnel  number.  For  n  =  5,  d/a  =  0.106  and  for  n  =  10,  d/a  =  0.051. 
The  curves  shown  in  Figs.  1  through  6  correspond  to  n  =  5. 

Before  comparing  finite-difference  solutions  with  fresnel-integral  solutions,  we  look  at  the  effect  of  apodization 
of  the  aperture  edge  for  the  fresnel-integral  solutions.  The  field  is  represented  on  a  square  grid  with  256  by  256 
points.  The  guard  ratio  is  rg  =  1.1.  That  is,  the  grid  width  is  1.1  times  the  width  of  the  initial  aperture.  The 
Phelps  algorithm  does  not  suffer  from  aliasing.  The  number  of  points  over  which  the  field  is  smoothed  at  the 
aperture  edge,  by  a  cosine  window,  varies  for  the  plots  in  Fig.  1.  The  number  of  points  is:  Curve  1,  nw  =  0; 
Curve  2,  nw  =  5;  Curve  3,  nw  =  10;  Curve  4,  nw  =  15.  The  first  three  curves  are  very  similar  to  each  other,  and 
the  fourth  is  noticeably  smoother.  The  outer  fresnel  zone  has  a  width  corresponding  to  25  points.  For  all  of  the 
fresnel-integral  results  shown  here,  a  grid  of  256  by  256  points  is  used,  and  rg  =  1.1  except  as  otherwise  noted. 


X 


Figure  1.  Calculated  intensity  for  the  propagation  of  a  beam  with  uniform  intensity,  that  has  passed  through  a  round 
aperture.  The  fresnel  number  is  5.  The  kernel  averaged  fresnel-integral  algorithm  was  used.  The  beam  was  smoothed 
with  a  cosine  window  over  nw  points.  For  Curve  1,  nw  =  0;  for  Curve  2,  nw  =  5;  for  Curve  3,  nw  =  10;  for  Curve  4, 
nw  =  15. 


Curves  shown  in  Fig.  2  compare  intensity  plots  for  a  fresnel-integral  solution  (solid)  and  for  the  finite- 
difference  solution  (dashed)  for  nw  =  5.  For  the  finite-difference  soultion,  rg  =  2.0  and  the  grid  was  466  by 
466  points,  so  that  the  same  number  of  points  were  inside  the  initial  aperture  for  both  calculations.  The  plots 
are  very  similar,  with  minor  differences  in  the  small-scale  structure.  The  physical  problem  is  not  exactly  the 
same  for  the  two  solutions.  The  fresnel  integral  corresponds  to  a  field  propagating  through  unobstructed  space, 
whereas  the  finite-difference  includes  a  cylindrical  absorbing  wall  with  twice  the  initial  diameter  of  the  beam. 

Because  the  grid  doesn’t  match  the  symmetry  of  the  optical  beam,  it  is  worth  checking  how  well  the  cylindrical 
symmetry  of  the  initial  aperture  is  preserved  after  propagation  using  the  finite-difference  algorithm.  For  the  finite- 
difference  calculation  of  Fig.  2,  the  intensity  was  evaluated  by  interpolation,  along  radial  lines  at  angles,  with 
the  a:-axis,  of  0,  7r / 12,  7t/6,  and  7r/4.  The  four  curves,  shown  in  Fig.  3,  are  almost  identical. 

Curves  shown  in  Fig.  4  are  intensity  plots  for  a  fresnel-integral  solution  (solid),  with  nw  =  10,  and  for  a 
finite-difference  solution  (dashed).  For  the  finite-difference  soultion,  the  grid  was  232  by  232  points,  half  as  many 
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Figure  2.  Intensity  calculated  with  the  fresnel-integral  algorithm,  solid  line,  and  the  intensity  calculated  with  the  finite- 
difference  algorithm,  dashed  line.  The  aperture  contained  the  same  number  of  transverse  points  for  each  of  the  two 
algorithms,  and  the  fresnel  number  was  5. 
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Figure  3.  Field  calculated  with  the  finite-difference  algorithm,  for  a  fresnel  number  of  5.  Traces  of  the  intensity,  evaluated 
by  interpolation,  are  shown  for  angles,  with  the  x-axis,  of  0  Rad  (Curve  1),  7t/12  Rad,  7r/6  Rad,  and  7t/4  Rad  (Curve  4). 


in  each  dimension  as  for  Fig.  2,  The  initial  aperture  was  smoothed  over  nw  =  5  points,  so  that  d/a  would  be  the 
same  for  the  two  curves.  The  plots  are  not  as  close  to  each  other  as  for  Fig.  2. 

Figure  5  shows  plots  of  the  intensity  for  a  grid  with  232  by  232  points.  As  for  Fig.  3,  the  intensity  was 
evaluated  by  interpolation,  along  radial  lines  at  angles  ,with  the  x-axis,  of  0,  7r/12,  7 r/6,  and  7r/4.  Unlike  Fig.  3, 
these  plots  are  noticeably  different. 

In  Fig.  6,  calculations  with  different  numbers  of  grid  points  along  2  are  compared.  As  A2  becomes  larger 
than  a  certain  value,  we  expect  significant  error  to  occur.  We  represent  A 2  in  terms  of  the  dimensionless  variable 
a.  Values  of  a  are  1  ,2  ,4,  and  8  for  curves  1  through  4,  respectively.  The  curve  for  a  =  8  is  quite  different  from 
the  others. 

Plots  for  a  fresnel  number  of  10  are  shown  in  Fig.  7.  The  pairs  of  curves  show  the  intensity  calculated 
with  the  fresnel-integral  algorithm,  solid  line,  and  the  intensity  calculated  with  the  finite-difference  algorithm, 
dashed  line.  The  finite-difference  (FD)  solution  on  the  bottom  has  a  resolution  of  466  by  466  grid  points.  The 
corresponding  curve  of  the  upper  pair  has  a  resolution  of  932  by  932  points.  For  the  bottom  pair,  nw  =  5,  and 
for  the  upper  FD  solution,  nw  =  10.  The  curve  with  higher  resolution  is  in  substantially  better  agreement  with 
the  fresnel-integral  solution. 
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Figure  4.  Intensity  calculated  with  the  fresnel-integral  algorithm,  solid  line,  and  the  intensity  calculated  with  the  finite- 
difference  algorithm,  dashed  line.  The  transverse  resolution  for  the  finite-difference  calculation  was  half  as  great  as  for 
Fig.  2. 


Figure  5.  Field  calculated  with  the  finite-difference  algorithm,  for  a  fresnel  number  of  5.  Traces  of  the  intensity,  evaluated 
by  interpolation,  are  shown  for  angles,  with  the  x-axis,  of  0  Rad  (Curve  1),  n/12  Rad,  n/6  Rad,  and  7r/4  Rad  (Curve  4). 
The  transverse  resolution  for  the  finite-difference  calculation  was  half  as  great  as  for  Fig.  3. 


4.  COMPARISON  OF  UNSTABLE  RESONATOR  MODES 

As  final  examples,  we  calculate  the  empty  cavity  modes  for  unstable  resonators.  The  examples  are  for  resonators 
with  magnification  of  2.  Equivalent  collimated  propagation6, 7  was  applied  for  the  entire  round  trip  (from  the 
outcoupling  plane  back  to  the  same  plane)  for  the  finite-difference  calculation,  and  interpolation  was  used  to 
magnify  the  beam.  The  fresnel-integral  (FI)  solution  was  for  a  confocal  unstable  resonator  with  a  back  mirror 
that  filled  the  grid  whose  width  was  1.2  times  the  geometric  beam  width.  Equivalent  collimated  propagation 
was  applied  for  each  of  the  two  propagation  steps,  and  interpolation  was  also  used  for  this  case,  to  magnify  the 
beam.  For  both  examples:  nw  =  5,  the  FI  grid  had  256  by  256  points,  and  the  FD  grid  had  466  by  466  points. 

Plots  of  the  intensity  for  an  unstable  resonator  with  an  equivalent  fresnel  number  of  2.3  are  shown  in  Fig.  8. 
The  solid  line  corresponds  to  the  fresnel-integral  algorithm,  and  the  dashed  line  resulted  from  the  finite-difference 
algorithm.  For  this  case,  agreement  is  rather  good,  although  the  curves  differ  significantly  near  the  resonator 
axis.  The  magnitudes  calculated  for  the  eigenvalue  were  0.616  and  0.625,  respectively. 

Plots  of  the  intensity  for  an  unstable  resonator  with  an  equivalent  fresnel  number  of  4.3  are  shown  in  Fig.  9. 
The  solid  line  corresponds  to  the  fresnel-integral  algorithm,  and  the  dashed  line  resulted  from  the  finite-difference 
algorithm.  The  FD  grid  had  466  by  466  points.  Agreement  is  not  good.  The  transverse  resolution  is  apparently 
inadequate  to  resolve  the  field  variation  for  this  case. 
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Figure  6.  Field  calculated  with  the  finite- difference  algorithm,  for  a  fresnel  number  of  5.  The  solid  curves  are  for  a  =  0.49. 
The  dashed  traces  are  for  varous  values  of  a:  a  —  1.0,  Curve  1;  a  =  2.0,  Curve2;  a  =  4.0,  Curve  3;  and  a  =  8.0,  Curve  4. 
The  transverse  grid  had  466  by  466  points. 


Figure  7.  Intensity  calculated  with  the  fresnel-integral  algorithm,  solid  line,  and  the  intensity  calculated  with  the  finite- 
difference  algorithm,  dashed  line.  The  finite-difference  solution  on  the  bottom  has  a  resolution  of  466  by  466  grid  points. 
The  corresponding  curve  of  the  upper  pair  has  a  resolution  of  932  by  932  points.  The  fresnel  number  was  10. 


5.  CONCLUDING  REMARKS 

The  effects  of  aperture  diffraction  can  be  treated  adequately  using  the  alternating-direction  implicit  finite- 
difference  algorithm,  although  transverse  resolution  requirements  are  rather  high.  For  the  cases  considered 
here,  the  approximate  requirement  for  resolution  in  the  z-direction  is  given  by  a  <  4.  Unstable  resonator  modes 
can  also  be  obtained  using  the  finite-difference  algorithm,  if  the  problem  is  adequately  resolved  by  the  number 
of  transverse  grid  points. 

Fresnel-integral  solutions  are  only  valid  if  <f>  is  a  small  angle, 

(j)  =  tan-1  ^  (10) 

where  l  is  the  propagation  distance.  If  this  condition  is  not  met,  finite-difference  solutions  can  not  be  validated 
by  comparison  with  them.  For  propagation  distances  much  shorter  than  this,  geometrical  optics,  with  some  spill 
over  at  the  beam  edges,  provides  a  reasonable  approximation. 

It  is  useful  to  understand  the  parameter  values  needed  to  treat  aperture  diffraction.  The  simulation  of 
lasers  usually  involves  compromise  for  the  sake  of  practicality,  however.  As  was  observed  by  Anan’ev,8  in  lasers 
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Figure  8.  Intensity  for  unstable  resonator  calculated  with  the  fresnel-integral  algorithm,  solid  line,  and  the  intensity 
calculated  with  the  finite-difference  algorithm,  dashed  line.  The  equivalent  fresnel  number  was  2.3,  and  the  resolution 
for  the  finite-difference  solution  was  466  by  466  points.  The  intensity  is  shown  for  the  field  incident  on  the  outcoupling 
mirror.  The  edge  of  the  outcoupling  mirror  is  at  x=l. 


Figure  9.  Intensity  for  unstable  resonator  calculated  with  the  fresnel-integral  algorithm,  solid  line,  and  the  intensity 
calculated  with  the  finite-difference  algorithm,  dashed  line.  The  intensity  is  shown  for  the  held  incident  on  the  outcoupling 
mirror.  The  equivalent  fresnel  number  was  4.3,  and  the  resolution  for  the  finite-difference  solution  was  466  by  466  points. 


with  large  outcoupling,  which  is  usually  the  case  if  unstable  resonators  are  used,  the  beam  shape  is  determined 
mostly  by  the  gain  distribution  rather  than  by  diffraction.  Precise  treatment  of  aperture  diffraction  may  play  a 
secondary  role,  whereas  the  wave-optics  treatment  of  strong  aberrations  may  be  quite  significant  in  determining 
the  laser  beam  quality.  It  is  probably  acceptable  to  use  a  finite-difference  propagation  algorithm,  with  a  number 
of  transverse  grid  points  that  would  be  inadequate  for  the  calculation  of  bare  cavity  modes,  to  predict  the 
power  and  beam  quality  of  a  laser  with  an  unstable  resonator.  Often,  the  gain  medium  is  separated  by  a 
considerable  distance  from  one  or  both  resonator  end  mirrors.  Then  it  would  be  efficient  to  combine  fresnel- 
integral  propagation  for  empty  regions  with  finite-difference  propagation  in  regions  with  strong  aberrations  or 
nonlinear  behavior. 
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